home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / LIBIP / FLTRGAUS.C < prev    next >
C/C++ Source or Header  |  1999-09-11  |  4KB  |  118 lines
  1. /* 
  2.  * fltrgaus.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. #include <stdlib.h>
  10. #include <malloc.h>
  11. #include <math.h>
  12.  
  13. #define ROOTLOG2 0.832554611    /* for Gaussian calc of std. dev. */
  14.  
  15. /*
  16.  * fltrgaus()
  17.  *   DESCRIPTION:
  18.  *     fltrgaus performs Gaussian filtering in frequency domain
  19.  *     with specified low & high pass bounds
  20.  *   ARGUMENTS:
  21.  *     imgReal(float **) pointer to real image array
  22.  *     imgImag(float **) pointer to imaginary image array
  23.  *     nRow, nCol(long) number of rows and columns for real and img arrays
  24.  *     passType(short)  passType=0 low pass
  25.  *                              passType=1 high pass
  26.  *                              passType=2 band pass
  27.  *                              passType=3 stop pass
  28.  *     f1, f2(double) low, high pass frequency bounds
  29.  *   RETURN VALUE:
  30.  *     none
  31.  */
  32.  
  33. void
  34. fltrgaus (imgReal, imgImag, nRow, nCol, passType, f1, f2)
  35.      float **imgReal, **imgImag;  /* real,complex image arrays */
  36.      long nRow, nCol;           /* size of image */
  37.      short passType;            /* lowpass=0, high=1, band=2, stop=3 */
  38.      double f1, f2;             /* low,high pass frequency bounds */
  39. {
  40.   register long x, y, nRowD2, nColD2,  /* half of row/column */
  41.     nRowM1, nColM1;             /* row/col minus 1 */
  42.   double bigger;                /* passband fraction fct of bigger axis */
  43.   double *fltr;                 /* filter vector */
  44.   double stdDevL, stdDevH;      /* std. dev.s of low and high pass fltrs */
  45.   long maxRad;                  /* max. sampling radius */
  46.   long dist;                    /* euclidean distance to point */
  47.   double tempL, tempH;
  48.   long i;
  49.  
  50. /* allocate filter vector -- radius from (u,v)=(0,0) to corners of freq plot */
  51.   nRowD2 = nRow / 2;
  52.   nColD2 = nCol / 2;
  53.   bigger = (double) (nRowD2 > nColD2) ? nRowD2 : nColD2;
  54.   maxRad = (long) sqrt ((double) (nRowD2 * nRowD2 + nColD2 * nColD2)) + 1;
  55.   fltr = (double *) calloc (maxRad, sizeof (double));
  56.  
  57.  
  58. /* determine filter vector coefficients for chosen type of filter */
  59.   switch (passType) {
  60.   case 0:                      /* low-pass */
  61.     stdDevL = f1 / ROOTLOG2 * bigger;
  62.     for (i = 0; i < maxRad; i++)
  63.       fltr[i] = exp (-(double) i * (double) i / (2.0 * stdDevL * stdDevL));
  64.     break;
  65.   case 1:                      /* high-pass */
  66.     stdDevH = f1 / ROOTLOG2 * bigger;
  67.     for (i = 0; i < maxRad; i++)
  68.       fltr[i] = 1.0 -
  69.         exp (-(double) i * (double) i / (2.0 * stdDevH * stdDevH));
  70.     break;
  71.   case 2:                      /* band-pass */
  72.     stdDevH = f1 / ROOTLOG2 * bigger;
  73.     stdDevL = f2 / ROOTLOG2 * bigger;
  74.     for (i = 0; i < maxRad; i++) {
  75.       tempL = exp (-(double) i * (double) i / (2.0 * stdDevL * stdDevL));
  76.       tempH = 1.0 -
  77.         exp (-(double) i * (double) i / (2.0 * stdDevH * stdDevH));
  78.       fltr[i] = tempL * tempH;
  79.     }
  80.     break;
  81.   case 3:                      /* and band-stop */
  82.     stdDevL = f1 / ROOTLOG2 * bigger;
  83.     stdDevH = f2 / ROOTLOG2 * bigger;
  84.     for (i = 0; i < maxRad; i++) {
  85.       tempL = exp (-(double) i * (double) i / (2.0 * stdDevL * stdDevL));
  86.       tempH = 1.0 -
  87.         exp (-(double) i * (double) i / (2.0 * stdDevH * stdDevH));
  88.       fltr[i] = tempL + tempH;
  89.     }
  90.     break;
  91.   default:
  92.     break;
  93.   }
  94.  
  95. /* perform filtering */
  96.   nRowM1 = nRow - 1;
  97.   nColM1 = nCol - 1;
  98.   for (y = 0; y < nRowD2; y++) {
  99.     for (x = 0; x < nColD2; x++) {
  100.       dist = (long) (sqrt ((double) (y * y + x * x)) + 0.5);
  101.  
  102.  
  103.       imgReal[y][x] *= (float) fltr[dist];
  104.       imgReal[y][nColM1 - x] *= (float) fltr[dist];
  105.       imgReal[nRowM1 - y][x] *= (float) fltr[dist];
  106.       imgReal[nRowM1 - y][nColM1 - x] *= (float) fltr[dist];
  107.  
  108.       imgImag[y][x] *= (float) fltr[dist];
  109.       imgImag[y][nColM1 - x] *= (float) fltr[dist];
  110.       imgImag[nRowM1 - y][x] *= (float) fltr[dist];
  111.       imgImag[nRowM1 - y][nColM1 - x] *= (float) fltr[dist];
  112.     }
  113.   }
  114.  
  115.  
  116.   return;
  117. }
  118.